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Abstract 



We study the ground state of a system of Bose hard-spheres trapped in 
an isotropic harmonic potential to investigate the effect of the interatomic 
correlations and the accuracy of the Gross-Pitaevskii equation. We compare 
a local density approximation, based on the energy functional derived from 
the low density expansion of the energy of the uniform hard sphere gas, and 
a correlated wave function approach which explicitly introduces the correla- 
tions induced by the potential. Both higher order terms in the low density 
expansion, beyond Gross-Pitaevskii, and explicit dynamical correlations have 
effects of the order of percent when the number of trapped particles becomes 
similar to that attained in recent experiments (N~ 10 7 ). 
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The recent discovery of Bose-Einstein condensation (BEC) of magnetically trapped alkali 
atoms has generated a huge amount of theoretical investigations. A review of the present 
situation can be found in refs.01. Present experimental conditions are such that the atomic 
gas is very dilute, i.e., the average distance among the atoms is much larger than the 
range of the interaction. As a consequence, the physics should be dominated by two body 
collisions, generally well described in terms of the s-wave scattering length. The case of 
a positive scattering length is equivalent to consider a very dilute system of hard spheres 
(HS), whose diameter coincides with the scattering length itself. So, the Gross-Pitaevskii 
(GP) theory for weakly interacting bosons seems the logical tool to study these systems and 
most of the present days theoretical work is founded on (or has its starting point in) this 
theory!. However, in very recent experiments the number of trapped atoms has spectacularly 
increased reaching N values of the order 10 6 and 10 7 atomsi. Therefore, it seems logical to 
ask for a deeper study of the effect of the interatomic correlations and of the accuracy of 
the GP scheme in this new scenario. 

In order to address these questions, we study the ground state of a system of Bose hard 
spheres trapped by a harmonic potential. More precisely, we consider hard spheres with 
a diameter of 52. 9A, corresponding to the s-wave triplet-spin scattering length of 87 Rb, in 
an isotropic harmonic trap characterized by an angular frequency uj/2tt = 77.78Hz. We 
also examine the large N atomic sodium case of ref.Q. We use and compare two methods: 
(i) a local density approximation (LDA) based on an energy functional derived by the low- 
density expansion of the energy of an uniform hard sphere gas; (ii) a correlated basis function 
(CBF) approach which explicitely takes into account the dynamical correlations induced by 
the potential and which is not, in principle, limited to purely repulsive interactions. 

LDA theory. In the case of an uniform hard sphere Bose gas, a perturbation theory in 

■ 3 where n is the density of the system, leads to the following 
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The energy functional associated with the GP theory is simply obtained in LDA by 
keeping only the first term in the expansion ([!]): 
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where the wavef unction ip is normalized to the total number of atoms. 

By performing a functional variation of Egp[i/j] one finds the Euler-Lagrange equation, 
known as the GP equation 
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where \x is the chemical potential. This equation has the form of a nonlinear stationary 
Schrodinger equation, and it has been solved for several types of traps using different nu- 
merical methodsiHi. 
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A logical step further, in the spirit of LDA, is to include into the energy functional 
the next terms of the correlation energy (JJ). It is convenient to simplify the notation 
by expressing lengths and energies in harmonic oscillator (HO) units. The HO length is 
a ho — (h/muj) 1 / 2 , and the spatial coordinates, the energy and the wave function are rescaled 
as r = auo^i E = huEi and ip(r) = (N / a HO ) 1 ^ 2 ^i(f) , where ifti(f) is normalized to unity 

By taking into account the next terms of the correlation energy, LDA provides a modified 
GP (MGP) energy functional 
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where a = a/ano and fii is the chemical potential in HO units. The MGP equation contains 
extra nonlinar terms in ipi. 

CBF approach. CBF theorv is a powerful tool to study strongly interacting many-body 
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systems (for a review, see ref.EI). In particular, it was used to study the HS homogenous 
Fermi gas in order to ascertain the accuracy of low density type expansions^. 
For N interacting bosons, at T=0 temperature, described by the hamiltonian 

H = ~ E V 2 + E *Ur 4 ) + E Vinj), (6) 
zm i i i<j 

where V ex t(^i) is the confining potential and V(ry) is the interatomic potential, the CBF 
ground-state wave function is 

i/>cbf(1, -,N) = F(l, .., N)1> MF (1, .., N). (7) 

where F(l, .., N) is a many-body correlation operator applied to the mean field w.f. ipMF- 
The advantage of using a CBF basis lies in the fact that non perturbative effects, as the 
short range repulsion for the hard spheres, may be directly inserted into the correlation. 
The simplest correlation operator has the Jastrow formE! 

F(l,..,JV) = n/j(r«), (8) 

i<j 

where the Jastrow correlation function, fj(r), depends on the interparticle distance 
only. fj(r) is variationally determined by minimizing the expectation value of Ecbf = 
{^cbf\H\%Ijcbf) I ' (4>cbf\iPcbf) ■ Eqbf may be evaluated either by MonteCarlo techniques 
or by cluster expansions of the needed one- and two-body densities, pi(vi) and p 2 (ri,i"2)- 
The energy per particle can be written as E CBF jN = T p + V e + V corr , where 
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and the correlation energy, V corr , is 
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where p\ is normalized to N. 

p 2 may be calculated by cluster expansion in powers of the dynamical correlation, h(r) = 
fj(r) — 1, and the integral Hypernetted Chain (HNC) equations^ may be used to sum infinite 
classes of cluster diagrams. 

Given the low-density of the trapped bosons system, it looks plausible to start by a lowest 
order (LO) cluster expansion. In this approximation p 2 LO \rx,r 2 ) = pi(ri)pi(r 2 )/j(ri 2 ) and 

Vg$ = ~ J dv x J dr 2 pi LO \r u r 2 )V JF (r 12 ) (12) 

where Vjf{t) = V{r) + (h 2 /m)[W In fj(r)] 2 is the Jackson-Feenberg potential. 

The minimization of E^bf with respect to p\ leads to the LO-correlated Hartree 
(CH/LO) equation 



p 1 1 /2 (r)=pp\ /2 (r) (13) 



with s — |r — rjj. 

The optimal choice for the Jastrow factor would be the one satisfying the Euler equa- 
tion 5EcBF/5fj(r) = 0. Otherwise, parametrized functional forms may be chosen whose 
parameters are found through the minimization process. We adopt here the correlation 
function minimizing the lowest order energy of a homogenous Bose gas with a healing con- 
dition at a distance d (taken as variational parameter). For the HS case, fj(r < a) = and 
fj(r > a) = u(r)/r where u(r) is solution of the Schrodinger-like equation — u" = K 2 u. 
fj(r) has the formB 

dmn[K(r-a)] 

JA ) rsm[K(d-a)} { V 

where the healing conditions, fj(r>d) = l and fj(r — d) — 0, are fixed by the relation: 
cot[K(d - a)] = (Kd)- 1 . 

Results. We begin by briefly discussing the boson HS homogeneous case. It was shown 
in ref.M that the lowest order cluster energy of an infinite system of bosonic hard spheres, 
correlated by the Jastrow factor flH[), asymptotically tends to the first term of the LD 
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expansion (LDO) in (0) when the healing distance goes to infinity. We have numerically 
checked this behavior. The situation changes if all the many-body cluster terms are included 
via HNC summation, since the computed energy shows a clear minimum in d. Some results 
are given in Fig.l for various x- values. The energies have been multiplied by lO 3 ^ 1 ) at 
x = lCT 6 ^ 5 ' -4 ), respectively. The Figure shows the energy estimates computed by retaining 
different LD expansion terms (LD1 and LD2 correspond to adding the second and third 
terms in (ID)), the HNC energies and, at the highest x-value, the exact energy, evaluated 
by a diffusion MonteCarlo methodEl. The quality of the HNC results is excellent as they 
practically coincide with the exact ones at all the densities (the low density value at low x 
and the DMC one at large x). The convergence of the LD expansion at large x appears to 
be rather poor, pointing to the relevance of successive contributions. 

The GP, MGP and CH/LO equations for ipi have been solved by the steepest descent 
method^ for the isotropic harmonic oscillator trap previously described. An initial trial 
state is projected onto the minimum of the functional by propagating it in imaginary time. 
In practice, one chooses an arbitrary time step At and iterates the equation 

^(rx, t + At) « ^(n, t) - AtiJVi(ri, t) (15) 

by normalizing ipi to 1 at each iteration. When the number of particles becomes very 
large the time step, which governs the rate of convergence, should be taken accordingly 
small. The number of iterations substantially increases and it is convenient to start the 
convergence process from a reasonable wave function (for the GP and MGP equations we 
start from the ground state harmonic oscillator function which minimizes the GP energy 
functional, while, in the CH/LO case, we start from the GP solution). 

We find that the LO approximation in the finite system shows a behavior similar to the 
infinite case. In particular, the solution of the CHE/LO equation provides an energy (in 
reduced unities) very close to the GP one when the healing distance becomes very large. 
In fact, for N = 10 5 and a = 4.33 x 10 -3 , corresponding to the 87 Rb scattering length, 
we obtain e[ lo) /N = 12.57, 12.28 and 12.11 at d = 10, 12 and 15, respectively, while 
E[ GP ^ /N = 12.10. This indicates that many-body effects can be recovered only by a HNC 
treatment also for a finite number of atoms. However, the solution of the fully correlated 
HNC Hartree equation, even if feasible in principle, is very cumbersome, computationally 
time consuming and numerically delicate. For this reason we have decided to estimate these 
effects by adopting a LDA-type approach to V corr . We approximate V corr ~ V£R A where 
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^imc(Pi) i s the HNC homogeneous gas energy per particle at density p\. 

The minimization of the energy gives the HNC correlated Hartree equation (CH/HNC) 



^i(r) = fJ.itpi(r) (17) 



where we have again introduced the scaled unities and x = pitr = Nd ■ 

The calculations have been performed for the 87 Rb scattering length. The scaled energies 
per particle and the root mean square radii are reported in Table I for particle numbers from 
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10 3 to 10 7 . The Table also shows the results obtained by neglecting the kinetic energy term 
in the GP equation. This approach, lousely called Thomas Fermi (TF) approximation, has 
been discussed in the literature and allows for deriving simple analytical expressions!!. The 
differences between this TF approach and a rigorous one have been recently discussed0ll! 
for spatially inhomogeneous Bose condensates. LDA has been used0 to estimate corrections 
to the GP energy within the TF approximation and retaining only the first correction in 
(|I|). The second correction is negative and partially cancels the first one. For instance, the 
cancellations go from ~ 15% for N = 10 4 to ~ 40% at iV = 10 6 if we just take the central 
densities, whereas the final energy is reduced by ~ 15% at N = 10 6 and it is practically 
unaffected by the second correction at lower iV-values. 

As expected, the TF results are close to the GP ones when N becomes large. The 
differences between GP and MGP increase with the number of particles and are of the order 
of 4% for the chemical potential and 2.5% for the energy at iV = 10 7 . The higher order 
terms in the low density expansion always have a repulsive effect. The same behavior is 
shown by the HNC results, which, however, are less repulsive than MGP at the large N 
values. 

We notice that if one uses the GP solution to perturbatively estimate the MGP energy, 
then the correction is negative (at N = 10 7 , AEi = —4.54). The non linear character of 
the MGP equation is responsible for this discrepancy. 

The density profile (normalized to unity) for N = 10 7 particles is given in Fig. 2. For 
this large number of particles the TF and GP densities are close, whereas the more repulsive 
MGP and HNC solutions lower the central density, expanding the density distribution and 
providing a larger radius, as shown in Table I . 

We have also considered a system of N = 1.5 X 10 Na atoms (a = 27.5 A) in a spherical 
trap having a frequency of 230 Hz. These conditions roughly correspond to those of the 
experiment described in ref.i. The results are shown in the last row of the Table and in 
Fig. 2. The effects of the correlations are similar to those found in the large N Rb cases. 
The energy increases by ~ 1% and the r.m.s. radius by ~ 0.7% respect to GP. The HNC 
central density is slightly reduced. 

In conclusion, we find that both higher order terms in the low density expansion (beyond 
the GP approach and evaluated in LDA) and explicit dynamical correlations (induced by the 
strong repulsion) have a not negligible effect when the number of trapped particles becomes 
large. So, they must be properly considered. 

In this respect, CBF theory may play a preminent role. In addition, it allows for a fully 
microscopic investigation for any type of potential. This is particularly interesting in view of 
its application to atoms having negative scattering lengths, as 7 Li. In this case the potential 
exhibits an actractive part and the GP equation has a metastable solution only if N does not 
exceed some critical value N c . A recent studyEl, which makes use of an effective interaction, 
has shown that a new branch of Bose condensate may appear at higher densities. CBF may 
provide further insights into this problem having access to the full structure of the potential 
itself. 

The introduction of a HNC energy functional derived by the homogeneous HS system 
provides a quick and probably reliable way of embodying correlation effects into the treat- 
ment of the trapped atoms. We particular appeal of this approach the clear possibilty 
of extending it to non spherical traps, as in real experiments. Moreover, as already stated, 
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potentials other than the simple hard sphere one may be easily handled. 
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TABLES 



TABLE I. Ground state energies per particle and chemical potentials of iV 87 Rb atoms in an 
isotropic trap (uj/2-k = 77.78 Hz). The last row refers to the Na case (u/2ir = 230 Hz). Energies 
in units of Tiu) and radii in HO lengths. 
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2.66 


3.04 


3.06 


3.04 


1.90 


2.43 


2.43 


2.43 


1.65 


1.66 


1.66 


10 4 


6.67 


6.87 


6.92 


6.89 


6.87 


5.04 


5.08 


5.04 


2.44 


2.45 


2.44 


10 5 


16.75 


16.85 


17.07 


16.94 


11.96 


12.10 


12.25 


12.20 


3.80 


3.84 


3.83 


10 6 


42.07 


42.12 


42.97 


42.53 


30.05 


30.12 


30.66 


30.48 


6.01 


6.10 


6.06 


10 7 


105.68 


105.70 


108.75 


107.20 


75.49 


75.52 


77.48 


76.85 


9.52 


9.74 


9.64 


1.5xl0 7 


91.07 


91.10 


92.41 


91.67 


65.05 


65.09 


65.92 


65.66 


8.84 


8.92 


8.90 
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FIGURES 



FIG. 1. Energy per particle (in units of h 2 /2ma 2 ) for hogeneous hard spheres in function of 
x. See text. 

FIG. 2. Density profiles for N = 10 7 Rb atoms and for N = 1.5 x 10 7 Na atoms in different 
approaches (dotted line=GP, dashed line=MGP, solid line=HNC). Densities are normalized to 
unity and distances are in units of HO lengths. 
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